Phase transition in site-diluted Josephson-junction arrays: A numerical study 
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Intriguing effects produced by random percolative disorder in two-dimensional Josephson-junction 
arrays are studied by means of large-scale numerical simulations. Using dynamic scaling analysis, we 
evaluate critical temperatures and critical exponents in high accuracy. With the introduction of site- 
diluted disorder, the Kosterlitz-Thouless phase transition is eliminated and evolves into continuous 
phase transition with a power-law divergent correlation length. Moreover, genuine depinning transi- 
tion and the related creep motion are studied, distinct types of creep motion for different disordered 
systems are observed. Our results not only are in good agreement with the recent experimental 
findings, but also shed some light on the relevant phase transitions. 
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I. INTRODUCTION 



Understanding the critical behavior of Josephson- 
junction arrays (JJA's) with various disorders is always 
a challenging question and has been intensely studied 
in recent years [l]]-[l(j. However, the properties of dif- 
ferent phases and various phase transitions are not well 
understood. Josephson-junction arrays gives an excellent 
realization to both two-dimensional (2D) XY model and 
granular High-T c superconductors . As is well known 
that the pure JJA's undergoes celebrated Kosterlitz- 
Thouless (KT) phase transition [ijj from the supercon- 
ducting state to the normal state, this transition is driven 
by the unbinding of thermally created topological defects. 
When the disorder is introduced, the interplays between 
the periodic pinning potential caused by the discreteness 
of the arrays, the repulsive vortex- vortex interaction and 
the effects produced by the disorder provide a rich phys- 
ical picture. 

In diluted JJA's, islands are randomly removed from 
the square lattice. Since it is a representative model for 
realizing the irregular JJA's systems, how the percola- 
tion influences the physical properties of JJA's has at- 
tracted considerable attention p], 0, 0, EL @- Harris et 
al introduced random percolative disorder into Nb-Au- 
Nb proximity-coupled junctions, the current-voltage (7- 
V) characteristics were measured and the results demon- 
strated that the only difference of the phase transition 
compared with that in ideal JJA's system is the decrease 
of critical temperature, while the phase transition still 
belong to the KT-type with the disorder strength span- 
ning from p — 0.7 to p — 1.0 (here 1 — p is the fraction of 
diluted sites) However, a recent experimental study 
by Yun et al showed that the KT-type phase transition in 
unfrustrated JJA's was eliminated due to the introduc- 
tion of site-diluted disorder [9( . Therefore, the existence 
of the KT-type phase transition in site-diluted JJA's re- 
mains a topic of controversy. 

On the other hand, much attention has been paid 
to investigate the zero-temperature depinning transition 
and the related low-temperature creep motion both the- 



oretically [lj, L!J, Bjj and numerically [1J, [lj, in a 
large variety of physical problems, such as charge density 
waves [l3j], random- field Ising model fl6j | . and flux lines 
in type-II superconductors [TtI [l8[ . Since the non-linear 
dynamic response is a striking aspect, there is increasing 
interest in these systems , espec ially in the flux lines of 
type-II superconductors [13, [3 • In a recent numerical 
study on the three-dimensional glass states of flux lines, 
Arrhenius creep motion and non-Arrhenius creep motion 
were observed with strong collective pinning and weak 
collective pinning, respectively [l7| . 

In this work, we numerically investigate the finite- 
temperature phase transition in site-diluted JJA's at dif- 
ferent percolative disorders, the zero-temperature depin- 
ning transition and the low-temperature creep motion are 
also considered. The outline of this paper is as follows. 
Section[TT]describes the model and the numerical method 
briefly. In section HTT1 we present the simulation results, 
analyzing them by means of scaling analysis. Section UVl 
gives a short summary of the main conclusions. 



II. MODEL AND SIMULATION METHOD 

JJA's can be described by the 2D XY model on a sim- 
ple square lattice, the Hamiltonian of which is [l^, [2(| 
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where the sum is over all nearest neighboring pairs on a 
2D square lattice, Jy denotes the strength of Josephson 
coupling between site i and site j, 4>i specifies the phase 
of the superconducting order parameter on site i, and 
Ay = (27r/$o) / A ■ d\ is the integral of magnetic vector 
potential from site i to site j with $o the flux quantum. 
The direct sum of Ay around an elementary plaquette is 
27r/, with / the magnetic flux penetrating each plaquette 
produced by the uniformly applied field, measured in unit 
of $o- In this paper, / = and / = 2/5 are in focus. 
The system sizes are selected as 128 x 128 for / = and 
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100 x 100 for / = 2/5, the finite size effects in these sizes 
are negligible. Diluted sites are randomly selected, then 
the nearest four bonds of which are removed from the 
lattice. The same random-number seed is used to choose 
the diluted sites, the percolative threshold concentration 
is about 0.592(1) for both systems [HI]. 

The resistivity-shunted-junction (RSJ) dynamics is in- 
corporated in the simulations, which can be described as 

dam 

— 2^ - <f>j ) = - — + J ex ,i - 2^ thj , (2) 

3 1 3 

where a is the normal conductivity, J ex ^ refers to the 
external current, rjij denotes the thermal noise current 
with < r)ij(t) >= and < 7?y (*)%(*') >= 2ak B TS(t - 
t). 

The fluctuating twist boundary condition is applied in 
the xy plane to maintain the current, thus the new phase 
angle 9i = <fii + ri ■ A (A = (A x , A y ) is the twist variable) 
is periodic in each direction. In this way, supercurrent 
between site i and site j is given by J"->j = Jijsin(0i — 
9j — Aij — nj ■ A) , and the dynamics of A a can be written 
as 



^ a = J2 [ Ji ^' + Vij] ~ I ^ (3) 

<i,j>a 

where a denotes the x or y direction, the voltage drop 
in a direction is V — —LA a . For convenience, units are 
taken as 2e = h = Jq = a = ks = 1 in the following. 
Above equations can be solved efficiently by a pseudo- 
spectral algorithm due to the periodicity of phase in all 
directions. The time stepping is done using a second- 
order Runge-Kutta scheme with At = 0.05. Our runs are 
typically (4 — 8) x 10 7 time steps and the latter half time 
steps are for the measurements. The detailed procedure 
in the simulations was described in Ref. [2(J H2] • In this 
work, a uniform external current I along x direction is 
fed into the system. 

Since RSJ simulations with direct numerical integra- 
tions of stochastic equations of motion are very time- 
consuming, it is practically difficult to perform any se- 
rious disorder averaging in the present rather large sys- 
tems. Our results are based on one realization of disor- 
der. For these very large samples, it is expected to exist 
a good self-averaging effect, which is confirmed by two 
additional simulations with different realizations of dis- 
order. This point is also supported by a recent study 
of JJA's by Um et al they confirmed that a well- 
converged disorder averaging for the measurement is not 
necessary, and well-converged data for large systems at 
a single disorder realization leads to a convincing result. 
In addition, simulations with different initial states are 
performed and the results are independent on the initial 
state we used. Actually, the hysteric phenomenon is usu- 
ally negligible in previous RSJ dynamical simulations on 



JJA's 0,1- For these reasons, the results from simula- 
tions with a unique initial state (random phases in this 
work) are accurate and then convincing. 

III. RESULTS AND DISCUSSION 
A. Finite temperature phase transition 

The I-V characteristics are measured at various disor- 
der strengths and temperatures. At each temperature, 
we try to probe the system at a current as low as possi- 
ble. To check the method used in this work, we investi- 
gate the I-V characteristics for / = 0,p = 1.0. As shown 
in Fig. 1(a), the slope of the I-V curve in log-log plot 
at the transition temperature T c w 0.894 is equal to 3, 
demonstrating that the I-V index jumps from 3 to 1, in 
consistent with the well-known fact that the pure JJA's 
experiences a KT-type phase transition at T c sw 0.894. 
Figs. 1(b) and (c) show the I-V traces at different per- 
colative disorders in unfrustrated JJA's, while Fig. 1(d) 
for / = 2/5, p = 0.65. It is clear that, at lower temper- 
atures, R — V/I tends to zero as the current decreases, 
which follows that there is a true superconducting phase 
with zero linear resistivity. 

It is crucial to use a powerful scaling method to an- 
alyze the I-V characteristics. In this paper, we adopt 
the Fisher-Fisher-Huse (FFH) dynamic scaling method, 
which provides an excellent approach to analyze the 
superconducting phase transition (23j. If the properly 
scaled I-V curves collapse onto two scaling curves above 
and bellow the transition temperature, a continuous 
superconducting phase transition is ensured. Such a 
method is widely used recently [f| [24| , the scaling form 
of which in 2D is 

v=ic z mio, (4) 

where V'+f— )( x ) is the scaling function above (below) T c , 
z is the dynamic exponent, £ is the correlation length, 
and V - I z+1 at T = T c . 

Assuming that the transition is continuous and char- 
acterized by the divergence of the characteristic length 
£ ~ \T - T c \~ v and time scale t ~ £ z , FFH dynamic 
scaling takes the following form 

(V/I)\T - T c \~ z » = yj ± (I\T - T c \~ v ). (5) 

On the other hand, a new scaling form is successfully 
adopted to certify a KT-type phase transition in JJA's 
by O 

(I/T)(I/V)V*=P ± (IZ/T), (6) 

note that the Eq. ([6]) can be obtained directly from the 
FFH dynamic scaling form after some simple algebra. 
The correlation length of KT-type phase transition above 
T c is well defined as £ ~ e ( c /l T ~ Tc D 1/2 and Eq. © reads 

(I/T)(I/V)^ Z = P+(Ie^ T - T ^ 1/2 /T). (7) 
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FIG. 1: I-V characteristics for different frustrations and per- 
colative disorders. The dash lines are drawn to show where 
the phase transition occurs, the slopes of which are equal to 
z + 1, z is the dynamic exponent. The transition temperature 
and dynamic exponent for (a) are well consistent with the 
well-known result, i. e., T c = 0.894, z = 2.0, for (b),(c),(d) 
are well consistent with those determined by FFH dynamic 
scaling analysis. Solid lines are just guide to eyes. 



As shown in Fig. 2, using T c = 0.24 ± 0.01, z = 
1.2 ± 0.02 and v = 1.0 ± 0.02, we get an excellent col- 
lapse for f = 0,p = 0.65 according to equation ([5]). In 
addition, all the low-temperature I-V curves can be fit- 
ted to V ~ Iexp(-(a/J)' 4 ) with (i = 0.9 ~ 1.1. These 
results certify a continuous superconducting phase with 
long-rang phase coherence. The critical temperature for 
such a strongly disordered system is very close to that in 
2D gauge glass model (T = 0.22) 

For / = 0,p = 0.86, firstly, we still adopt the scaling 
form in equation ([5]) to investigate the I-V character- 
istics. As displayed in Fig. 3, we get a good collapse 
for T < T c with T c = 0.58 ± 0.01, z = 2.0 ± 0.01 and 
v = 1.4 ± 0.02, demonstrating a superconducting phase 
with power-law divergent correlation for T < T c . Note 
that the collapse is bad for T > T c , indicating that the 
phase transition is not a completely non-KT-type one. 
Next, we use the scaling form in equation to ana- 
lyze the I-V data above the critical temperature. In- 
terestingly, using T c = 0.58 and z = 2.0 determined 
above, a good collapse for T > T c is achieved, which is 
shown in Fig. 4. That is to say, the I-V characteristics 
at T < T c are like those of a continuous phase transi- 
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FIG. 2: Dynamic scaling of I-V data at various temperatures 
according to equation ([5]) for / = 0, p = 0.65. 



tion with power-law divergent correlation length while 
at T > T c are like those of KT-type phase transition, 
which are well consistent with the recent experimental 
observations [?|. Therefore, by the present model, we 
recover the phenomena in experiments and give some 
insight into the phase transition. More information on 
the low-temperature phase calls for further equilibrium 
Monte Carlo simulations as in Ref. (27| . 

To make a comprehensive comparison with the exper- 
imental findings as in Ref. [9j, we also investigated the 
finite-temperature phase transition in frustrated JJA's 
(/ = 2/5) at a strong site-diluted disorder (p = 0.65). As 
shown in Fig. 5, a superconducting phase transition with 
power-law divergent correlation is clearly observed. As is 
well known, non-KT-type phase transition in frustrated 
systems is a natural result. However, it is intriguing to 
see that in unfrustrated systems, one may ask what our 
results really imply and what is the mechanism for it. It 
has been revealed that in the presence of a strong random 
pinning which is produced by random site dilutions, a 
breaking of ergodicity due to large energy barrier against 
vortex motion may allow enough vortices to experience a 
non-KT-type continuous transition [28[ . 

The systems considered in our work are site-diluted 
JJA's, which are not the same as bond-diluted JJA's in 
Ref. [all. The difference is, in bond-diluted systems the 
diluted bonds are randomly removed, while in the site- 
diluted systems, the diluted sites are randomly selected, 
then the nearest four bonds around the selected sites are 
removed. Although the JJA's in Ref. [H, 0| and the 
present work are diluted in different ways, it is interesting 
to note that some of the obtained exponents are very 
close, possibly due to the similar disorder effect produced. 



B. Depinning transition and creep motion 

Next, we pay attention to the zero-temperature depin- 
ning transition and the related low-temperature creep 
motion for the typical site-diluted JJA's systems men- 
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FIG. 3: Dynamic scaling of I-V data at various temperatures 
according to equation ([5]) for / = 0,p = 0.86, T < T c . Solid 
lines are just guide to eyes. 



FIG. 5: Dynamic scaling of I-V data at various temperatures 
according to equation ([5]) for / = 2/5, p — 0.65. 
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FIG. 4: Dynamic scaling of I-V data at various temperatures 
according to equation (|7| for / = 0,p = 0.86, T > T c . Solid 
lines are just guide to eyes. 



tioned above. Depinning can be described as a critical 
phenomenon with scaling law V ~ (I — I C )P, demonstrat- 
ing a transition from a pinned state below critical driving 
force I c to a sliding state above I c . The (7 — I c ).vs.V 
traces at T = for / = 0,p = 0.86; / = Q,p = 0.65 and 
/ = 2/5, p = 0.65 are displayed in Fig. 6, linear-fittings 
of Log(I — I c ).vs.LogV curves are also shown as solid 
lines. As for / = 0,p = 0.86, the depinning exponent (3 
is determined to be 2.62 ± 0.1 and the critical current I c 
is 0.302 ± 0.005, while for the cases f = Q,p= 0.65 and 
/ = 2/5, p = 0.65, the depinning exponents are evaluated 
to be 2.37 ±0.1 and 2.27 ±0.05 with the critical currents 
I c = 0.039 ± 0.001 and I c = 0.035 ± 0.002, respectively. 

When the temperature increases slightly, creep mo- 
tions can be observed. In the low-temperature regime, 
the I-V traces are rounded near the zero-temperature 
critical current due to thermal fluctuations. Fisher first 
suggested to map such a phenomenon for the ferromag- 
net in magnetic field where the second-order phase tran- 
sition occurs [5^1 . This mapping was then extended to 
the random-field Ising model [161 ] and the flux lines in 
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FIG. 6: {a)IV characteristics for / = 0,p = 0.86 with I c = 
0.302 ± 0.005, P = 2.62 ± 0.1. (bjlV characteristics for / = 
0,p = 0.65 with I c = 0.039 ± 0.001, = 2.37 ± 0.1. (c)IV 
characteristics for / = 2/5, p = 0.65 with I c = 0.035 ± 0.002, 
P = 2.27 ±0.05 
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FIG. 7: (a)LogV-LogT curves for / = 0,p = 0.86 in the 
vicinity of I c with I c = 0.302 ± 0.001, 1/8 = 1.688 ± 0.001. 
(b)LogV-LogT curves for f — 0,p = 0.65 in the vicinity of 
I c with I c = 0.03875 ± 0.0005, 1/5 = 2.24 ± 0.02. (c)LogV- 
LogT curves for / = 2/5, p = 0.65 in the vicinity of I c with 
Ic = 0.034 ± 0.001, 1/5 = 2.29 ± 0.01. 



type-II superconductors [l7J . For the flux lines in type-II 
superconductors, if the voltage is identified as the order 
parameter, the current and the temperature are taken as 
the inverse temperature and the field respectively, anal- 
ogous to the second-order phase transition in the ferro- 
magnet, the voltage, current and the temperature will 
satisfy the following scaling ansatz [13, [HI 



V(T, I) = T^ S S[{1 - Ic/^T- 1 ^ 5 }. 



(8) 



The relation V(T, I = I c ) = 5 , (0)T 1 / 5 can be easily de- 
rived at / = I c , by which the critical current I c and the 
critical exponent 8 can be determined through the linear 
fitting of the LogT — LogV curve at I c . The LogT — LogV 
curves are plotted in Fig. 7(a) for / = 0,p = 0.86. 
We can observe that the critical current is between 0.3 
and 0.32. In order to locate the critical current pre- 
cisely, we calculate other values of voltage at current 
within (0.3,0.32) with a current step 0.01 by quadratic 
interpolation [261 ]. Deviation of the T-V curves from 
the power law is calculated as the square deviations 



SD = ^2[V(T) — y(T)] 2 between the temperature range 
we calculated, here the functions y(T) = C1T~ C2 are ob- 
tained by linear fitting of the LogT — LogV curves. The 
current at which the SD is minimum is defined as the 
critical current. The critical current is then determined 
to be 0.302 ± 0.001. Simultaneously, we obtain the expo- 
nent 1/8 = 1.688 ± 0.001 from the slope of LogT - logV 
curve at I c = 0.302. The similar method is applied to in- 
vestigate the cases / = 0,p = 0.65 and / = 2/5, p = 0.65. 
As shown in Figs. 7 (b) and (c), the critical current I c 
and critical exponent 1/8 for f = 0,p = 0.65 are deter- 
mined to be 0.03875 ± 0.0005, 2.24 ± 0.02 respectively, 
for / = 2/5, p = 0.65, the result is I c = 0.034 ± 0.001, 
1/8 = 2.29 ±0.01. 

We then draw the scaling plots according to Eq. 8. 
Using the one parameter tuning of (3, we get the best 
collapses of data to a single scaling curve with (3 = 
2.61 ± 0.02 and 2.28 ± 0.02 for / = 0,p = 0.86 and 
/ = 0,p = 0.65 in the regime I < I c , respectively, which 
are shown in Figs. 8(a) and (b). For / = 0,p = 0.86, this 
curve can be fitted by S(x) — 0.0994exp(1.9a:), combined 
with the relation j3S — 1.55, suggesting a non-Arrhcnius 
creep motion. However, for the strongly site-diluted sys- 
tem with / = Q,p = 0.65, the scaling curve can be fitted 
by S(x) — 0.037exp(0.5a;), combined with the relation 
[38 w 1.0, indicative of an Arrhenius creep motion. Inter- 
estingly, as displayed in Fig. 8(c) for / = 2/5, p = 0.65, 
the exponent (3 is fitted to be 2.30 ± 0.02, which yields 
(38 « 1.0. The scaling curve in the regime I < I c can be 
fitted by S{x) = 0.105ea;p(0.25a;). These two combined 
facts suggest an Arrhenius creep motion in this case. 

It is worthwhile to note that both the finite- 
temperature phase transition and the creep motion for 
strongly disordered JJA's (p = 0.65) with and with- 
out frustration are very similar. The TV curves in 
low temperature for all three cases can be described by 
V oc T^-^expiAil - I C /I)/T 0S ), which is just one of the 
main characteristics of glass phases [13, Hij]. While the 
TV traces for KT-type phases can be fitted to V oc I a . 
Therefore, we have provided another evidence for the ex- 
istence of non-KT-type phases in the low-temperature 
regime for these three cases (/ = 0,p = 0.86; f = 0,p = 
Q.65;/ = 2/5,p = 0.65). 



IV. SUMMARY 

To explore the properties of various phase transitions 
in site-diluted JJA's, we have performed large scale simu- 
lations at two typical percolative strengths p = 0.86 and 
p = 0.65 as in a recent experimental work [9J]. The RSJ 
dynamics was incorporated in our work, from which we 
measured the TV characteristics at different tempera- 
tures. The critical temperature of the finite-temperature 
phase transition was found to decrease as the diluted 
sites increase. For / = 0,p = 0.86, the phase transition 
is the combination of a KT-type transition and a con- 
tinuous transition with power-law divergent correlation 
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FIG. 8: (a) Scaling plot for / = 0,p = 0.86 with I c = 0.302, 
1/5 = 1.688 and (35 = 1.55. (b) Scaling plot for / = 0,p = 
0.65 with Ic = 0.03875, 1/5 = 2.24 and (35 » 1.0. (c)Scaling 
plot for / = 2/5, p = 0.65 with I c = 0.034, 1/5 = 2.29 and 
/3<5 « 1.0. 

length. At strong percolative disorder (p = 0.65), the 
KT-type phase transition in pure JJA's is changed into a 
completely non-KT-type phase transition, moreover, the 
finite-temperature phase transition for frustrated JJA's 
is similar to that in unfrustrated JJA's. All the obtained 
dynamic exponents z = a — 1, with a the I-V index 
at the critical temperature, and all the static exponents 
fall in the range of v — (1.0,2.0) usually observed at 
vortex-glass transitions experimentally. Following table 



summarizes the critical temperatures at different frustra- 
tions and disorder strengths. 



TABLE I: Summary of T c 





f=0 


f=2/5 


p=0.95 


0.85(2) 


0.16(2) 


p=0.86 


0.58(1) 


0.13(1) 


p=0.7 


0.27(2) 


0.12(1) 


p=0.65 


0.24(1) 


0.14(1) 



In a recent experiment, Yun et al Q suggested a non- 
KT-type phase transition in unfrustrated JJA's with site- 
diluted disorder for the first time, however the nature of 
these phase transitions and various phases is still in an 
intensive debate. Our results not only recover the re- 
cent experimental findings but also shed some light 
on the various phases. Non-KT-type finite-temperature 
phase transition in site-diluted JJA's was confirmed by 
the scaling analysis. The different divergent correlations 
at various disorder strengths were suggested, the crit- 
ical exponents were evaluated in high accuracy, which 
are crucial for understanding such a critical phenomenon. 
Furthermore, the results in this paper are not only useful 
for understanding the site-diluted systems, but also use- 
ful for understanding the whole class of disordered JJA's. 
For instance, the combination of two different phase tran- 
sitions may exist in other disordered JJA's systems. 

In addition, the zero-temperature depinning transition 
and the low-temperature creep motion are also touched. 
It is demonstrated by the scaling analysis that the creep 
law for f = 0,p= 0.86 is non-Arrhenius type while those 
for / = Q,p = 0.65 and / = 2/5, p — 0.65 belong to 
the Arrhenius type. The evidence of non-KT-type phase 
transition can also be provided by this scaling analysis. It 
is interesting to note that the non-Arrhenius type creep 
law for weak disorder (f = 0,p — 0.86) is similar to 
that in three-dimensional flux lines with a weak collective 
pinning [T7| . The product of the two exponents 1.55 is 
also very close to 3/2 determined in Ref.[l7j|- For / = 
0,p = 0.65 and / = 2/5, p = 0.65, the observed Arrhenius 
type creep law is also similar to that in the glass states 
of flux lines with a strong collective pinning as in Ref. 

Future experimental work is needed to clarify this 
observation. 



V. ACKNOWLEDGEMENTS 

This work was supported by National Natural Sci- 
ence Foundation of China under Grant Nos. 10774128, 
PCSIRT (Grant No. IRT0754) in University in China, 
National Basic Research Program of China (Grant 
Nos. 2006CB601003 and 2009CB929104), and Zhejiang 
Provincial Natural Science Foundation under Grant No. 
Z7080203. 

' Corresponding author. Email:qhchen@zju. edu.cn 



7 



[1] D. C. Harris, S. T. Herbert, D. Stroud and J. C. Garland, 

Phys. Rev. Lett. 67, 3606 (1991). 
[2] E. Granato and D. Dominguez, Phys. Rev. B 56, 14671 

(1997). 

[3] M. Bcnakli and E. Granato, S. R. Shenoy and M. Gabay, 

Phys. Rev. B 57, 10314 (1998). 
[4] E . Granato and D. Dominguez, Phys. Rev. B 63, 094507 

(2001). 

[5] Y. J. Yun, I. C. Baek and M. Y. Choi, Phys. Rev. Lett. 

89, 037004 (2002). 
[6] E. Granto and D. Dominguez, Phys. Rev. B 71, 094521 

(2005). 

[7] J. S. Lim, M. Y. Choi, B. J. Kim and J. Choi, Phys. Rev. 

B 71, 100505R (2005). 
[8] J. Um, B. J. Kim, P. Minnhagen, M. Y. Choi and S. I. 

Lee, Phys. Rev. B 74, 094516 (2006). 
[9] Y. J. Yun, I. C. Baek and M. Y. Choi, Phys. Rev. Lett. 

97, 215701 (2006). 
[10] Y. J. Yun, I. C. Baek and M. Y. Choi, Europhys. Lett. 

76, 271 (2006). 

[11] C. J. Lobb, D. W. Abraham and M. Tinkham, Phys. 
Rev. B 27, 150 (1983); M. Prester, Phys. Rev. B 54, 606 
(1996). 

[12] J. M. Kosterlitz and D. J. Thouless, J.Phys.C 6, 1181 

(1973); J. M. Kosterlitz, J. Phys. C 7, 1046 (1974); V. 

L. Berezinskii, Sov. Phys. - JETP, 34,610 (1972); V. L. 

Berezinskii, Zh. Eksp. Teor. Fiz. 61,1144 (1973). 
[13] T. Nattermann, Phys. Rev. Lett. 64, 2454 (1990). 
[14] P. Chauve, T. Giamarchi and P. L. Doussal, Phys. Rev. 

B. 62, 6241 (2000). 
[15] M. Midler, D. A. Gorokhov and G. Blatter, Phys. Rev. 

B. 63, 184305 (2001). 



[16] L. Rosters, A. Hucht, S. Liibeck, U. Nowak and K. D. 

Usadel, Phys. Rev. E. 60, 5202 (1999). 
[17] M. B. Luo and X. Hu, Phys. Rev. Lett. 98, 267002 (2007). 
[18] P. Olsson, Phys. Rev. Lett. 98, 097001 (2007); Q. H. 

Chen, Phys. Rev. B 78, 104501 (2008). 
[19] P. Olsson and S. Teitel, Phys. Rev. Lett. 87, 137001 

(2001). 

[20] Q. H. Chen and X. Hu, Phys. Rev. Lett. 90, 117005 
(2003); Q. H. Chen and X. Hu, Phys. Rev. B 75, 064504 

(2007) . 

[21] T. Gebele, J. Phys. A: Math. Gen. 17, L51(1984); Y. 

Laroyer and E. Pommiers, Phys.Rev.B 50, 2795 (1994). 
[22] Q. H. Chen and L. H. Tang, Phys. Rev. Lett. 87, 067001 

(2001) ; L. H. Tang and Q. H. Chen, Phys. Rev. B 67, 
024508 (2003). 

[23] D. S. Fisher, M. P. A. Fisher and D. A. Huse, Phys. Rev. 

B 43, 130 (1991). 
[24] H. Yang, Y. Jia, L. Shan, Y. Z. Zhang, H. H. Wen, C. 

G. Zhuang, Z. K. Liu, Q. Li, Y. Cui and X. X. Xi, Phys. 

Rev. B 76, 134513 (2007). 
[25] J. Holzer, R. S. Newrock, C. J. Lobb, T. Aouaroun and 

S. T. Herbert, Phys. Rev. B 63, 184508 (2001). 
[26] Q. H. Chen, J. P. Lv, H. Liu, Phys. Rev. B 78, 054519 

(2008) . 

[27] H. G. Katzgraber, Phys. Rev. B 67, 180402R (2003); H. 
G. Katzgraber and A. P. Young, Phys. Rev. B 66 224507 

(2002) . 

[28] P. Holme and P. Olsson, Europhys. Lett. 60 439 (2002). 
[29] D. S. Fisher, Phys. Rev. Lett. 50, 1486 (1983); Phys. 
Rev. B 31, 1396 (1985). 



Phase transition in site-diluted Josephson-junction arrays: A numerical study 
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We numerically investigate the intriguing effects produced by random percolative disorder in 
two-dimensional Josephson-junction arrays. By dynamic scaling analysis, we evaluate critical tem- 
peratures and critical exponents with high accuracy. It is observed that, with the introduction 
of site-diluted disorder, the Kosterlitz-Thouless phase transition is eliminated and evolves into a 
continuous transition with power-law divergent correlation length. Moreover, genuine depinning 
transition and creep motion are studied, evidence for distinct creep motion types is provided. Our 
results not only are in good agreement with the recent experimental findings, but also shed some 
light on the relevant phase transitions. 
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I. INTRODUCTION 

Understanding the critical behavior of Josephson- 
junction arrays (JJA's) with various disorders is always 
a challenging issue and has been intensely studied in re- 
cent years [l|]-[l(j. However, the properties of different 
phases and various phase transitions are not well under- 
stood. Josephson-junction arrays gives an excellent re- 
alization to both two-dimensional (2D) XY model and 
granular High-T c superconductors [111 ]. As we know, the 
pure JJA's undergoes the celebrated Kosterlitz-Thouless 
(KT) phase transition from the superconducting state to 
the normal one, this transition is driven by the unbinding 
of thermally activated topological defects • When the 
disorder is introduced, the interplays among the repulsive 
vortex-vortex interaction, the periodic pinning potential 
caused by the discreteness of the arrays, and the defects 
produced by the disorder provide a rich physical picture. 

In site-diluted JJA's, the crosses around the randomly 
selected sites are removed from the square lattice. Since 
it is a representative model for realizing the irregular 
JJA's systems, how the percolation influences the physi- 
cal properties of JJA's has attracted considerable atten- 
tion [l], [3, [H HI Q ■ Harris et al introduced random per- 
colative disorder into Nb-Au-Nb proximity-coupled junc- 
tions, the current- voltage (I-V) characteristics were mea- 
sured and the results demonstrated that the only dif- 
ference of the phase transition compared with that in 
ideal JJA's system is the decrease of critical tempera- 
ture, while the transition type still belongs to the KT 
one with the disorder strength spanning from p = 0.7 
to p = 1.0 (here 1 — p is the fraction of diluted sites) 
However, in a recent experiment, Yun et al showed 
that the phase transition changes into a non-KT-type one 
when the disorder strength increases to a moderate value 
(p = 0.86) [9(. Therefore, the existence of the KT-type 
phase transition in site-diluted JJA's remains a topic of 
controversy, the nature of these phase transitions and the 
various phases is not clear. 

On the other hand, much effort has been devoted to 



the zero-temperature depinning transition (ZTDT) and 
the related low-temperature creep motion (LTCM) both 
theoretically [3 Hill and numerically [II G3, [l|] in a 
large variety of physical problems, such as charge density 
waves [l3[ , random- field Ising model (l6j , and flux lines in 
type-II superconductors Since the non- linear dy- 

namic response is a striking problem, there is increasing 
interest in its properties and characteristics, especially 
in the flux lines of type-II superconductors [13, LL8| • In 
a recent numerical study on the three-dimensional glass 
states of flux lines, Arrhcnius creep motion was observed 
at a strong collective pinning, while the non-Arrhenius 
creep motion was demonstrated at a weak collective pin- 
ning 

In this work, we numerically study the finite- 
temperature phase transition (FTPT) in site-diluted 
JJA's at different percolative disorder strengths, the 
ZTDT and the LTCM are also investigated. The outline 
of this paper is as follows. Sec. |TT] describes the model 
and the numerical method briefly. In Sec. IIII1 we present 
the main results, where some discussions are also made. 
Sec. IIVI gives a short summary of the main conclusions. 



II. MODEL AND SIMULATION METHOD 

JJA's can be described by the 2D XY model on a 
simple square lattice, the Hamiltonian reads [3, H3] 



H = - 



E 

<i,j> 



Jij cos( 



fa ~ A ij), 



(1) 



where the sum is over all nearest neighboring pairs on 
a 2D square lattice, Jij denotes the strength of Joseph- 
son coupling between site i and site j, fa specifies the 
phase of the superconducting order parameter on site i, 
and Aij = (27r/$o) J A • cil is the integral of magnetic 
vector potential from site i to site j , <&o denotes the flux 
quantum. The direct sum of around an elementary 
plaquette is 27r/, with / the magnetic flux penetrating 
each plaquette produced by the uniformly applied field, 
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which is measured in unit of <I>o- / = and / = 2/5 are 
the focuses of this paper. The system sizes are selected 
as 128 x 128 for / = and 100 x 100 for / = 2/5, where 
the finite size effects are negligible. We introduce the 
site-diluted disorder similar to the previous experiments 
P, Q . We first select the diluted sites randomly with the 
probability 1 — p, then remove the nearest four bonds 
around the selected sites from the lattice. The distribu- 
tions of the diluted sites are the same for all the samples 
considered. The percolative threshold concentration p c 
is about 0.592 (Hf. 

The resistivity-shunted-junction (RSJ) dynamics is in- 
corporated in the simulations, which can be described as 

turn 

<t h ^-^ ■ ■ dH ^-^ 

— 2J<& - <f>j) = -qT- + Jex,i - }^ na, (2) 
j 1 j 

where a is the normal conductivity, J ex ^ refers to the 
external current, r/ij denotes the thermal noise current 
with < T)ij(t) >= and < rj i3 (t)rj i: j(t r ) >= 2ak B T6(t - 
t). 

The fluctuating twist boundary condition is applied in 
the xy plane to maintain the current, thus the new phase 
angle 9i = <fii + ri ■ A (A = (A x , A y ) is the twist variable) 
is periodic in each direction. In this way, supercurrent 
between site i and site j is given by J'^j = Jijsin(9i — 
8j — Aij — rij • A) , and the dynamics of A Q can be written 
as 

<i,j>a 

where a denotes the x or y direction, the voltage drop 
in a direction is V = — LA a . For convenience, units are 
taken as 2e = h = Jq = a = ks = 1 in the following. 
Above equations can be solved efficiently by a pseudo- 
spectral algorithm due to the periodicity of phase in all 
directions. The time stepping is done using a second- 
order Runge-Kutta scheme with At = 0.05. Our runs are 
typically (4 — 8) x 10 7 time steps and the latter half time 
steps are for the measurements. The detailed procedure 
in the simulations was described in Ref . [2(J, H2| . In this 
work, a uniform external current I along x direction is 
fed into the system. 

Since RSJ simulations with direct numerical integra- 
tions of stochastic equations of motion are very time- 
consuming, it is practically difficult to perform any se- 
rious disorder averaging in the present rather large sys- 
tems. Our results are based on one realization of disor- 
der. For these very large samples, it is expected to exist 
a good self-averaging effect, which is confirmed by two 
additional simulations with different realizations of dis- 
order. This point is also supported by a recent study 
of JJA's by Um et al Q. In addition, simulations with 
different initial states are performed and the results are 
nearly the same. Actually, the hysteric phenomenon is 



usually negligible in previous RSJ dynamical simulations 
on JJA's [3, Q ■ For these reasons, the results from simu- 
lations with a unique initial state (random phases in this 
work) are accurate and convincing. 

III. RESULTS AND DISCUSSIONS 
A. Finite temperature phase transition 

The I-V characteristics are measured at different dis- 
order strengths and temperatures. At each temperature, 
we try to probe the system at a current as low as possi- 
ble. To check the method used in this work, we investi- 
gate the I-V characteristics for / = 0,p = 1.0. As shown 
in Fig. 1(a), the slope of the I-V curve in log-log plot 
at the transition temperature T c (« 0.894) is equal to 
3, demonstrating that the I-V index jumps from 3 to 1, 
consistent with the well-known fact that the pure JJA's 
experiences a KT-type phase transition at T c w 0.894. 
Figs. 1(b) and (c) show the I-V traces at different dis- 
order strengths in unfrustrated JJA's, while Fig. 1(d) is 
for / = 2/5, p = 0.65. It is clear that, at lower temper- 
atures, R — V/I tends to zero as the current decreases, 
which follows that there is a true superconducting phase 
with zero linear resistivity. 

It is crucial to use a powerful scaling method to an- 
alyze the I-V characteristics. In this paper, we adopt 
the Fisher-Fisher-Huse (FFH) dynamic scaling method, 
which provides an excellent approach to analyze the 
superconducting phase transition (23|. If the prop- 
erly scaled I-V curves collapse onto two scaling curves 
above and below the transition temperature, a continu- 
ous superconducting phase transition is ensured. Such a 
method is widely used @, [24[ , the scaling form of which 
in 2D is 

v = ic z MiO, (4) 

where ip + r_\(x) is the scaling function above (below) T c , 
z is the dynamic exponent, £ is the correlation length, 
and V - I z+1 at T = T c . 

Assuming that the transition is continuous and char- 
acterized by the divergence of the characteristic length 
f ~ |T - T c \~ v and time scale t ~ £ z , FFH dynamic 
scaling takes the following form 

(V/I)\T - T c \- Z » = i, ± (I\T - T c \~»). (5) 

On the other hand, to certify a KT-type phase tran- 
sition in JJA's, a new scaling form [25j is proposed as 
follows 

(I/T)(I/V)V*=P ± {IZ/T), (6) 

which can be derived directly from Eq. ^ after some 
simple algebra. The correlation length of KT-type phase 
transition above T c is well defined as £ ~ e ( c /l T - T d) 1/2 
and Eq. ^ is rewritten as 

(I/T^I/V) 1 ^ = P+(Ie^ T - T ^ 1/2 /T). (7) 
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FIG. 1: I-V characteristics for different frustrations and dis- 
order strengths. The dash lines are drawn to show where 
the phase transition occurs, the slopes of which are equal to 
z + 1, z is the dynamic exponent. The transition temperature 
and dynamic exponent for (a) are well consistent with the 
well-known result, i. e., T c = 0.894, z = 2.0, for (b),(c),(d) 
are well consistent with those determined by FFH dynamic 
scaling analysis. Solid lines are just guide to eyes. 



We perform the dynamic scaling analysis at a strong 
disorder (p — 0.65) in unfrustrated system (/ = 0). Us- 
ing T c = 0.24±0.01, z = 1.2±0.02 and v = 1.0±0.02, an 
excellent collapse is achieved according to Eq. ([5]) , which 
is shown in Fig. 2. In addition, all the low-temperature 
I-V curves can be fitted to V ~ iexp(— with 
/i = 0.9 ~ 1.1. These results certify a continuous super- 
conducting phase with long-rang phase coherence. The 
critical temperature for such a strongly disordered system 
is very close to that in 2D gauge glass model (T c = 0.22) 

For / = 0,p = 0.86, we first still adopt the scaling 
form in Eq. ([5]) to investigate the I-V characteristics. 
As displayed in Fig. 3, we get a good collapse for T <T C 
with T c = 0.58 ± 0.01, z = 2.0 ± 0.01 and u = lA± 0.02, 
demonstrating a superconducting phase with power-law 
divergent correlation for T < T c . Note that the collapse 
is poor for T > T c , implying that the phase transition 
is not a completely non-KT-type one. Next, we use the 
scaling form in Eq. to analyze the I-V data above 
T c . Interestingly, using T c = 0.58 and z = 2.0 determined 
above, a good collapse for T > T c is achieved, which is 
shown in Fig. 4. That is to say, the I-V characteristics 
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FIG. 2: Dynamic scaling of I-V data at various temperatures 
according to Eq. ([5]) for / = 0,p — 0.65. 



at T < T c are like those of a continuous phase transition 
with power-law divergent correlation length while at T > 
T c are like those of KT-type phase transition, which are 
well consistent with the recent experimental observations 

1- 

To make a comprehensive comparison with the experi- 
mental findings in Ref . 0] , we also investigate the FTPT 
in frustrated JJA's (/ = 2/5) at a strong site-diluted dis- 
order (p = 0.65). As shown in Fig. 5, a superconducting 
phase transition with power-law divergent correlation is 
clearly observed. 

As is well known, non-KT-type phase transition in 
frustrated systems is a natural result. However, it is in- 
triguing to see that in unfrustrated systems, one may 
ask what our results really imply and what is the mecha- 
nism. It has been revealed that in the presence of a strong 
random pinning which is produced by random site dilu- 
tions, a breaking of ergodicity due to large energy barrier 
against vortex motion may allow enough vortices to ex- 
perience a non-KT-type continuous transition (28j . 

Interestingly, we recover the phenomena in experi- 
ments by the present model and give some insights into 
the FTPT. More information on the low-temperature 
phase calls for further equilibrium Monte Carlo simu- 
lations as in Refs. [27| . Table I summarizes the crit- 
ical temperatures at different frustrations and disorder 
strengths. One can find that, the critical temperature 
in unfrustrated system decreases with increasing diluted 
sites. 



TABLE I: Summary of T c . 





f=0 


£=2/5 


p=0.95 


0.85(2) 


0.16(2) 


p=0.86 


0.58(1) 


0.13(1) 


p=0.7 


0.27(2) 


0.12(1) 


p=0.65 


0.24(1) 


0.14(1) 



The systems considered in our work are site-diluted 
JJA's, which are not the same as bond-diluted JJA's in 
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FIG. 3: Dynamic scaling of I-V data at various temperatures 
according to Eq. © for / = 0,p = 0.86, T < T c . Solid lines 
are just guide to eyes. 



FIG. 5: Dynamic scaling of I-V data at various temperatures 
according to Eq. ([5]) for / = 2/5, p = 0.65. 
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FIG. 4: Dynamic scaling of I-V data at various temperatures 
according to Eq. © for / = 0,p = 0.86, T > T c . Solid lines 
are just guide to eyes. 



Ref. H |H . In bond-diluted systems the diluted bonds 
are randomly removed, while in the site-diluted systems, 
the diluted sites are randomly selected, then the nearest 
four bonds around the selected sites are removed. Al- 
though the JJA's in Ref. [1, H[ and the present work are 
diluted in different ways, it is interesting to note that the 
obtained exponents in FTPT are very close, possibly due 
to the similar disorder effect produced. 



B. Depinning transition and creep motion 

Next, we turn to the ZTDT and the LTCM for the 
typical site-diluted JJA's systems mentioned above. De- 
pinning can be described as a critical phenomenon with 
scaling law V ~ (7 — ic) , demonstrating a transition 
from a pinned state below critical driving force I c to 
a sliding state above I c . The (I — I c ).vs.V traces at 
T = for / = 0,p = 0.86; / = 0,p = 0.65 and 
/ = 2/5, p = 0.65 are displayed in Fig. 6, linear-fittings 
of Log(I — I c ).vs.LogV curves are also shown as solid 
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FIG. 6: (a)/y characteristics for / = 0,p = 0.86 with I c = 
0.302 ± 0.005, P = 2.62 ± 0.1. (bjlV characteristics for / = 
0,p = 0.65 with I c = 0.039 ± 0.001, = 2.37 ± 0.1. (c)IV 
characteristics for / = 2/5, p = 0.65 with I c = 0.035 ± 0.002, 
f3 = 2.27 ±0.05. 
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satisfy the following scaling ansatz [17j, [26[ 

V(T, I) = T^ S S[{1 - IJI)^ 1 / 135 ]. 



(8) 
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FIG. 7: (a)LogT-LogV curves for / = 0,p = 0.86 around I c 
with I c = 0.302 ± 0.001, 1/8 = 1.688 ± 0.001. (b)LogT-LogV 
curves for / = 0,p = 0.65 around I c with I c = 0.03875 ± 
0.0005, 1/5 = 2.24 ±0.02. (c)LogT-LogV curves for / = 
2/5, p = 0.65 around I c with I c = 0.034 ±0.001, 1/8 = 2.29 ± 
0.01. 



lines. As for f = 0,p = 0.86, the depinning exponent (3 
is determined to be 2.62 ± 0.1 and the critical current I c 
is 0.302 ± 0.005, while for the cases f = 0,p= 0.65 and 
/ = 2/5, p = 0.65, the depinning exponents are evaluated 
to be 2.37 ±0.1 and 2.27 ±0.05 with the critical currents 
I c = 0.039 ± 0.001 and I c = 0.035 ± 0.002, respectively. 

When the temperature increases slightly, creep mo- 
tions can be observed. In the low-temperature regime, 
the I-V traces are rounded near the zero-temperature 
critical current due to thermal fluctuations. Fisher first 
suggested to map such a phenomenon for the ferromag- 
net in magnetic field where the second-order phase tran- 
sition occurs [29( | . This mapping was then extended to 
the random-held Ising model [161 ] and the flux lines in 
type-II superconductors [llj • For the flux lines in type-II 
superconductors, if the voltage is identified as the order 
parameter, the current and the temperature are taken as 
the inverse temperature and the field respectively, anal- 
ogous to the second-order phase transition in the ferro- 
magnet, the voltage, current and the temperature will 



Where S(x) is a scaling function. The relation V(T,I = 
I c ) = S f (0)r 1//(5 can be easily derived at I — I c , by which 
the critical current I c and the critical exponent 8 can be 
determined through the linear fitting of the LogT-Logy 
curve at I c . 

The LogT-LogV curves are plotted in Fig. 7(a) for 
/ = 0,p = 0.86. We can observe that the critical current 
is between 0.3 and 0.32. In order to locate the criti- 
cal current precisely, we calculate other values of voltage 
at current within (0.3,0.32) with a current step 0.01 by 
quadratic interpolation [26j]. Deviation of the T-V curves 
from the power law is calculated as the square deviations 
SD = ^2[V(T) — y(T)] 2 between the temperature range 
we calculated, here the functions y(T) = C1T~ C2 are 
obtained by linear fitting of the LogT-LogV/ curves. The 
current at which the SD is minimum is defined as the 
critical current. The critical current is then determined 
to be 0.302 ± 0.001. Simultaneously, we obtain the expo- 
nent 1/8 = 1.688 ± 0.001 from the slope of LogT-LogV 
curve at I c = 0.302. The similar method is applied to in- 
vestigate the cases / = 0,p = 0.65 and / = 2/5, p = 0.65. 
As shown in Figs. 7 (b) and (c), the critical current I c 
and critical exponent 1/8 for / = Q,p = 0.65 are deter- 
mined to be 0.03875 ± 0.0005, 2.24 ± 0.02 respectively, 
for / = 2/5, p = 0.65, the result is I c = 0.034 ± 0.001, 
1/8 = 2.29 ±0.01. 

We then draw the scaling plots according to Eq. ([5]). 
By one-parameter tuning (only (3), we get the best col- 
lapses of data in the regime I < I c with (3 = 2.61 ± 0.02 
and 2.28 ± 0.02 for / = 0,p = 0.86 and / = 0,p = 0.65 
respectively, which are shown in Figs. 8(a) and (b). For 
/ = Q,p = 0.86, this curve can be fitted by S(x) = 
0.0994exp(1.9a;), combined with the relation (38 — 1.55, 
suggesting a non-Arrhenius creep motion. However, for 
the strongly site-diluted system with f = 0,p = 0.65, 
the scaling curve can be fitted by S(x) = 0.037exp(0.5x), 
combined with the relation (38 « 1.0, indicative of an 
Arrhenius creep motion. Interestingly, as displayed in 
Fig. 8(c) for / = 2/5, p = 0.65, the exponent (3 is 
found to be 2.30 ± 0.02, which yields (38 w 1.0. The 
scaling curve in the regime I < I c can be fitted by 
S(x) = 0.105exp(0.25a;). These two combined facts sug- 
gest an Arrhenius creep motion in this case. 

It is worthwhile to note that both the FTPT and the 
LTCM for strongly disordered JJA's (p — 0.65) with and 
without frustration are very similar. The I-V curves 
in low temperature for all three cases can be described 
by V oc T 1 / 5 expL4(l - I c /I)/T^ s ], this is one of the 
main characteristics of glass phases [Tt], [Hj], while the 
I-V traces for KT-type phases can be fitted to V oc I a . 
Hence, by the scaling ansatze in Eq. ([5]), we have pro- 
vided another evidence for the existence of non-KT-type 
phases in the low-temperature regime for these three 
cases (/ = 0,p = 0.86;/ = Q,p = 0.65;/ = 2/5, p = 
0.65). 
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FIG. 8: (a) Scaling plot for / = 0,p = 0.86 with I c = 0.302, 
1/5 = 1.688 and (35 = 1.55. (b) Scaling plot for / = 0,p = 
0.65 with Ic = 0.03875, 1/5 = 2.24 and (35 » 1.0. (c)Scaling 
plot for / = 2/5, p = 0.65 with I c = 0.034, 1/5 = 2.29 and 
/3<5 « 1.0. 



IV. SUMMARY 

To explore the properties of various critical phenom- 
ena in site-diluted JJA's, we performed large-scale simu- 



lations for two typical percolative strengths p = 0.86 and 
p = 0.65 as in a recent experimental work 0]. We in- 
vestigated the FTPT, the ZTDT and the LTCM in these 
systems. The RSJ dynamics was applied in our work, 
from which we measured the I-V characteristics at dif- 
ferent temperatures. 

The results obtained in this work about FTPT are 
well consistent with the recent experimental findings in 
Ref. 0] and are inconsistent with the earlier experimen- 
tal study in Ref. [l[ , possibly due to the large noise in the 
measurement of voltage in Ref. [l[ (larger than 0.2nv), 
which was considerably reduced in the experiments by 
Yun et al 0] . The evidence for non-KT-type phase tran- 
sition was revealed by two different scaling ansatzes (Eq. 
(O and Eq. ©). Our results also shed some light on the 
various phases and the phase transitions where the differ- 
ent divergent correlations at various disorder strengths 
were suggested, and the critical exponents were evalu- 
ated. Furthermore, the results in this paper are useful 
for understanding not only the site-diluted systems, but 
also the whole class of disordered JJA's, for instance, the 
combination of two different phase transitions may exist 
in other disordered JJA's systems. 

In addition, the ZTDT and the LTCM were also 
touched. It was demonstrated by the scaling analysis 
that the creep law for / = Q,p = 0.86 is non-Arrhcnius 
type while those for / = 0,p = 0.65 and / = 2/5, p = 0.65 
belong to the Arrhenius type. It is interesting to note 
that the non- Arrhenius type creep law for / = 0,p = 0.86 
is similar to that in three-dimensional flux lines with a 
weak collective pinning [l7| • The product of the two ex- 
ponents 1.55 is also very close to 3/2 determined in Ref. 
|l7|]- For / = 0,p = 0.65 and / = 2/5, p = 0.65, the ob- 
served Arrhenius type creep law is also similar to that in 
the glass states of flux lines with a strong collective pin- 
ning as in Ref. [i"^ . Future experimental work is needed 
to clarify this observation. 
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We numerically investigate the intriguing effects produced by random percolative disorder in 
two-dimensional Josephson-junction arrays. By dynamic scaling analysis, we evaluate critical tem- 
peratures and critical exponents with high accuracy. It is observed that, with the introduction 
of site-diluted disorder, the Kosterlitz-Thouless phase transition is eliminated and evolves into a 
continuous transition with power-law divergent correlation length. Moreover, genuine depinning 
transition and creep motion are studied, evidence for distinct creep motion types is provided. Our 
results not only are in good agreement with the recent experimental findings, but also shed some 
light on the relevant phase transitions. 
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I. INTRODUCTION 

Understanding the critical behavior of Josephson- 
junction arrays (JJA's) with various disorders is always 
a challenging issue and has been intensely studied in re- 
cent years |l|-[lG|. However, the properties of different 
phases and various phase transitions are not well under- 
stood. Josephson-junction arrays gives an excellent re- 
alization to both two-dimensional (2D) XY model and 
granular High-T c superconductors As we know, the 
pure JJA's undergoes the celebrated Kosterlitz-Thouless 
(KT) phase transition from the superconducting state to 
the normal one, this transition is driven by the unbinding 
of thermally activated topological defects [u| • When the 
disorder is introduced, the interplays among the repulsive 
vortex-vortex interaction, the periodic pinning potential 
caused by the discreteness of the arrays, and the defects 
produced by the disorder provide a rich physical picture. 

In site-diluted JJA's, the crosses around the randomly 
selected sites are removed from the square lattice. Since 
it is a representative model for realizing the irregular 
JJA's systems, how the percolation influences the physi- 
cal properties of JJA's has attracted considerable atten- 
tion [3, [3, & EL 111 ■ Harris et al introduced random per- 
colative disorder into Nb-Au-Nb proximity-coupled junc- 
tions, the current- voltage (/-V) characteristics were mea- 
sured and the results demonstrated that the only dif- 
ference of the phase transition compared with that in 
ideal JJA's system is the decrease of critical tempera- 
ture, while the transition type still belongs to the KT 
one with the disorder strength spanning from p = 0.7 
to p = 1.0 (here 1 — p is the fraction of diluted sites) 
However, in a recent experiment, Yun et al showed 
that the phase transition changes into a non-KT-type one 
when the disorder strength increases to a moderate value 
(p = 0.86) @. Therefore, the existence of the KT-type 
phase transition in site-diluted JJA's remains a topic of 
controversy, the nature of these phase transitions and the 
various phases is not clear. 

On the other hand, much effort has been devoted to 
the zero-temperature depinning transition (ZTDT) and 



the related low-temperature creep motion (LTCM) both 
theoretically [HlllHl and numerically [II, El El in a 
large variety of physical problems, such as charge density 
waves [l3j , random- field Ising model [l6j , and flux lines in 
type-II superconductors [HEIj]- Since the non- linear dy- 
namic response is a striking problem, there is increasing 
interest in its properties and characteristics, especially 
in the flux lines of type-II superconductors [171 Il8j. In 
a recent numerical study on the three-dimensional glass 
states of flux lines, Arrhenius creep motion was observed 
at a strong collective pinning, while the non-Arrhenius 
creep motion was demonstrated at a weak collective pin- 
ning [H]. 

In this work, we numerically study the finite- 
temperature phase transition (FTPT) in site-diluted 
JJA's at different percolative disorder strengths, the 
ZTDT and the LTCM are also investigated. The outline 
of this paper is as follows. Sec. |TT] describes the model 
and the numerical method briefly. In Sec. IIII1 we present 
the main results, where some discussions are also made. 
Sec. IIVI gives a short summary of the main conclusions. 



II. MODEL AND SIMULATION METHOD 

JJA's can be described by the 2D XY model on a 
simple square lattice, the Hamiltonian reads [l!| [20| 



H = - 



E 

<i,j> 



Jij cos( 



0j — A-ij ) > 



(1) 



where the sum is over all nearest neighboring pairs on 
a 2D square lattice, Jij denotes the strength of Joseph- 
son coupling between site i and site j, fa specifies the 
phase of the superconducting order parameter on site i, 
and Aij = (27r/<I>o) / A • dl is the integral of magnetic 
vector potential from site i to site j , <&o denotes the flux 
quantum. The direct sum of Aij around an elementary 
plaquette is 27r/, with / the magnetic flux penetrating 
each plaquette produced by the uniformly applied field, 
which is measured in unit of $o- / — and f — 2/5 are 
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the focuses of this paper. The system sizes are selected 
as 128 x 128 for / = and 100 x 100 for / = 2/5, where 
the finite size effects are negligible. We introduce the 
site-diluted disorder similar to the previous experiments 
[H> @ • We first select the diluted sites randomly with the 
probability 1 — p, then remove the nearest four bonds 
around the selected sites from the lattice. The distribu- 
tions of the diluted sites are the same for all the samples 
considered. The_percolative threshold concentration p c 
is about 0.592 j^j. 

The resistivity-shuntcd-junction (RSJ) dynamics is in- 
corporated in the simulations, which can be described as 

turn 

ah sr-^ , • • . dH -r-^ . . 

IfelJ&i - = -qT + Jc " _ l^^ih ( 2 ) 

j 1 j 

where a is the normal conductivity, J ex ^ refers to the 
external current, rjij denotes the thermal noise current 
with < r]ij(t) >= and < %(*)%(*') >= 2ak B T5(t - 

t). 

The fluctuating twist boundary condition is applied in 
the xy plane to maintain the current, thus the new phase 
angle 9i = <fii + ri ■ A (A = (A x , A y ) is the twist variable) 
is periodic in each direction. In this way, supercurrent 
between site i and site j is given by Jf-^j = Jijsin(9i — 
8j — Aij — nj • A) , and the dynamics of A Q can be written 
as 



<i,j>a 

where a denotes the x or y direction, the voltage drop 
in a direction is V = — LA a . For convenience, units are 
taken as 2e = h = Jq = a = ks = 1 in the following. 
Above equations can be solved efficiently by a pseudo- 
spectral algorithm due to the periodicity of phase in all 
directions. The time stepping is done using a second- 
order Runge-Kutta scheme with At = 0.05. Our runs are 
typically (4 — 8) x 10 7 time steps and the latter half time 
steps are for the measurements. The detailed procedure 
in the simulations was described in Ref . [2(J, H2| ■ In this 
work, a uniform external current / along x direction is 
fed into the system. 

Since RSJ simulations with direct numerical integra- 
tions of stochastic equations of motion are very time- 
consuming, it is practically difficult to perform any se- 
rious disorder averaging in the present rather large sys- 
tems. Our results are based on one realization of disor- 
der. For these very large samples, it is expected to exist 
a good self-averaging effect, which is confirmed by two 
additional simulations with different realizations of dis- 
order. This point is also supported by a recent study 
of JJA's by Um et al In addition, simulations with 
different initial states are performed and the results are 
nearly the same. Actually, the hysteric phenomenon is 
usually negligible in previous RSJ dynamical simulations 



on JJA's 0,1- For these reasons, the results from simu- 
lations with a unique initial state (random phases in this 
work) are accurate and convincing. 

III. RESULTS AND DISCUSSIONS 
A. Finite temperature phase transition 

The I-V characteristics are measured at different dis- 
order strengths and temperatures. At each temperature, 
we try to probe the system at a current as low as possi- 
ble. To check the method used in this work, we investi- 
gate the I-V characteristics for / = 0,p = 1.0. As shown 
in Fig. 1(a), the slope of the I-V curve in log-log plot 
at the transition temperature T c (« 0.894) is equal to 
3, demonstrating that the I-V index jumps from 3 to 1, 
consistent with the well-known fact that the pure JJA's 
experiences a KT-type phase transition at T c sw 0.894. 
Figs. 1(b) and (c) show the I-V traces at different dis- 
order strengths in unfrustrated JJA's, while Fig. 1(d) is 
for / = 2/5, p = 0.65. It is clear that, at lower temper- 
atures, R — V/I tends to zero as the current decreases, 
which follows that there is a true superconducting phase 
with zero linear resistivity. 

It is crucial to use a powerful scaling method to an- 
alyze the I-V characteristics. In this paper, we adopt 
the Fisher-Fisher-Huse (FFH) dynamic scaling method, 
which provides an excellent approach to analyze the 
superconducting phase transition (23|. If the prop- 
erly scaled I-V curves collapse onto two scaling curves 
above and below the transition temperature, a continu- 
ous superconducting phase transition is ensured. Such a 
method is widely used @, [24| , the scaling form of which 
in 2D is 

v = ic z i>±(iO, (4) 

where V'+f— )( x ) is the scaling function above (below) T c , 
z is the dynamic exponent, £ is the correlation length, 
and V - I z+1 at T = T c . 

Assuming that the transition is continuous and char- 
acterized by the divergence of the characteristic length 
£ ~ \T - T c \- V and time scale t - £ z , FFH dynamic 
scaling takes the following form 

(V/I)\T - T c \- Z » = yj ± (I\T - T c \-'). (5) 

On the other hand, to certify a KT-type phase tran- 
sition in JJA's, a new scaling form [25j is proposed as 
follows 

(I/T)(I/V)^ =P ± m/T), (6) 

which can be derived directly from Eq. ([4]) after some 
simple algebra. The correlation length of KT-type phase 
transition above T c is well defined as £ ~ e ( c /l T - T d) 1 
and Eq. ([6]) is rewritten as 

(I/T^I/V) 1 ^ = P+(Ie^ T - T ^ 1/2 /T). (7) 
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FIG. 1: I-V characteristics for different frustrations and dis- 
order strengths. The dash lines are drawn to show where 
the phase transition occurs, the slopes of which are equal to 
z + 1, z is the dynamic exponent. The transition temperature 
and dynamic exponent for (a) are well consistent with the 
well-known result, i. e., T c = 0.894, z = 2.0, for (b),(c),(d) 
are well consistent with those determined by FFH dynamic 
scaling analysis. Solid lines are just guide to eyes. 



We perform the dynamic scaling analysis at a strong 
disorder (p — 0.65) in unfrustrated system (/ = 0). Us- 
ing T c = 0.24±0.01, z = 1.2±0.02 and v = 1.0±0.02, an 
excellent collapse is achieved according to Eq. ([5]) , which 
is shown in Fig. 2. In addition, all the low-temperature 
I-V curves can be fitted to V ~ iexp(— with 
/i = 0.9 ~ 1.1. These results certify a continuous super- 
conducting phase with long-rang phase coherence. The 
critical temperature for such a strongly disordered system 
is very close to that in 2D gauge glass model (T c = 0.22) 

For / = 0,p = 0.86, we first still adopt the scaling 
form in Eq. ([5]) to investigate the I-V characteristics. 
As displayed in Fig. 3, we get a good collapse for T <T C 
with T c = 0.58 ± 0.01, z = 2.0 ± 0.01 and u = lA± 0.02, 
demonstrating a superconducting phase with power-law 
divergent correlation for T < T c . Note that the collapse 
is poor for T > T c , implying that the phase transition 
is not a completely non-KT-type one. Next, we use the 
scaling form in Eq. to analyze the I-V data above 
T c . Interestingly, using T c = 0.58 and z = 2.0 determined 
above, a good collapse for T > T c is achieved, which is 
shown in Fig. 4. That is to say, the I-V characteristics 
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FIG. 2: Dynamic scaling of I-V data at various temperatures 
according to Eq. ([5]) for / = 0,p — 0.65. 



at T < T c are like those of a continuous phase transition 
with power-law divergent correlation length while at T > 
T c are like those of KT-type phase transition, which are 
well consistent with the recent experimental observations 

1- 

To make a comprehensive comparison with the experi- 
mental findings in Ref . 0] , we also investigate the FTPT 
in frustrated JJA's (/ = 2/5) at a strong site-diluted dis- 
order (p = 0.65). As shown in Fig. 5, a superconducting 
phase transition with power-law divergent correlation is 
clearly observed. 

As is well known, non-KT-type phase transition in 
frustrated systems is a natural result. However, it is in- 
triguing to see that in unfrustrated systems, one may 
ask what our results really imply and what is the mecha- 
nism. It has been revealed that in the presence of a strong 
random pinning which is produced by random site dilu- 
tions, a breaking of ergodicity due to large energy barrier 
against vortex motion may allow enough vortices to ex- 
perience a non-KT-type continuous transition (28j . 

Interestingly, we recover the phenomena in experi- 
ments by the present model and give some insights into 
the FTPT. More information on the low-temperature 
phase calls for further equilibrium Monte Carlo simu- 
lations as in Refs. [27| . Table I summarizes the crit- 
ical temperatures at different frustrations and disorder 
strengths. One can find that, the critical temperature 
in unfrustrated system decreases with increasing diluted 
sites. 



TABLE I: Summary of T c . 





f=0 


£=2/5 


p=0.95 


0.85(2) 


0.16(2) 


p=0.86 


0.58(1) 


0.13(1) 


p=0.7 


0.27(2) 


0.12(1) 


p=0.65 


0.24(1) 


0.14(1) 



The systems considered in our work are site-diluted 
JJA's, which are not the same as bond-diluted JJA's in 
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FIG. 3: Dynamic scaling of I-V data at various temperatures 
according to Eq. © for / = 0,p = 0.86, T < T c . Solid lines 
are just guide to eyes. 



FIG. 5: Dynamic scaling of I-V data at various temperatures 
according to Eq. ([5]) for / = 2/5, p = 0.65. 
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FIG. 4: Dynamic scaling of I-V data at various temperatures 
according to Eq. © for / = 0,p = 0.86, T > T c . Solid lines 
are just guide to eyes. 



Ref. H |H . In bond-diluted systems the diluted bonds 
are randomly removed, while in the site-diluted systems, 
the diluted sites are randomly selected, then the nearest 
four bonds around the selected sites are removed. Al- 
though the JJA's in Ref. [1, H[ and the present work are 
diluted in different ways, it is interesting to note that the 
obtained exponents in FTPT are very close, possibly due 
to the similar disorder effect produced. 



B. Depinning transition and creep motion 

Next, we turn to the ZTDT and the LTCM for the 
typical site-diluted JJA's systems mentioned above. De- 
pinning can be described as a critical phenomenon with 
scaling law V ~ (7 — ic) , demonstrating a transition 
from a pinned state below critical driving force I c to 
a sliding state above I c . The (I — I c ).vs.V traces at 
T = for / = 0,p = 0.86; / = 0,p = 0.65 and 
/ = 2/5, p = 0.65 are displayed in Fig. 6, linear-fittings 
of Log(I — I c ).vs.LogV curves are also shown as solid 
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FIG. 6: (a)/y characteristics for / = 0,p = 0.86 with I c = 
0.302 ± 0.005, P = 2.62 ± 0.1. (bjlV characteristics for / = 
0,p = 0.65 with I c = 0.039 ± 0.001, = 2.37 ± 0.1. (c)IV 
characteristics for / = 2/5, p = 0.65 with I c = 0.035 ± 0.002, 
f3 = 2.27 ±0.05. 
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satisfy the following scaling ansatz [17j, [26[ 

V(T, I) = T^ S S[{1 - IJI)^ 1 / 135 ]. 



(8) 
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FIG. 7: (a)LogT-LogV curves for / = 0,p = 0.86 around I c 
with I c = 0.302 ± 0.001, 1/8 = 1.688 ± 0.001. (b)LogT-LogV 
curves for / = 0,p = 0.65 around I c with I c = 0.03875 ± 
0.0005, 1/5 = 2.24 ±0.02. (c)LogT-LogV curves for / = 
2/5, p = 0.65 around I c with I c = 0.034 ±0.001, 1/8 = 2.29 ± 
0.01. 



lines. As for f = 0,p = 0.86, the depinning exponent (3 
is determined to be 2.62 ± 0.1 and the critical current I c 
is 0.302 ± 0.005, while for the cases f = 0,p= 0.65 and 
/ = 2/5, p = 0.65, the depinning exponents are evaluated 
to be 2.37 ±0.1 and 2.27 ±0.05 with the critical currents 
I c = 0.039 ± 0.001 and I c = 0.035 ± 0.002, respectively. 

When the temperature increases slightly, creep mo- 
tions can be observed. In the low-temperature regime, 
the I-V traces are rounded near the zero-temperature 
critical current due to thermal fluctuations. Fisher first 
suggested to map such a phenomenon for the ferromag- 
net in magnetic field where the second-order phase tran- 
sition occurs [29( | . This mapping was then extended to 
the random-held Ising model [161 ] and the flux lines in 
type-II superconductors [llj • For the flux lines in type-II 
superconductors, if the voltage is identified as the order 
parameter, the current and the temperature are taken as 
the inverse temperature and the field respectively, anal- 
ogous to the second-order phase transition in the ferro- 
magnet, the voltage, current and the temperature will 



Where S(x) is a scaling function. The relation V(T,I = 
I c ) = S f (0)r 1//(5 can be easily derived at I — I c , by which 
the critical current I c and the critical exponent 8 can be 
determined through the linear fitting of the LogT-Logy 
curve at I c . 

The LogT-LogV curves are plotted in Fig. 7(a) for 
/ = 0,p = 0.86. We can observe that the critical current 
is between 0.3 and 0.32. In order to locate the criti- 
cal current precisely, we calculate other values of voltage 
at current within (0.3,0.32) with a current step 0.01 by 
quadratic interpolation [26j]. Deviation of the T-V curves 
from the power law is calculated as the square deviations 
SD = ^2[V(T) — y(T)] 2 between the temperature range 
we calculated, here the functions y(T) = C1T~ C2 are 
obtained by linear fitting of the LogT-LogV/ curves. The 
current at which the SD is minimum is defined as the 
critical current. The critical current is then determined 
to be 0.302 ± 0.001. Simultaneously, we obtain the expo- 
nent 1/8 = 1.688 ± 0.001 from the slope of LogT-LogV 
curve at I c = 0.302. The similar method is applied to in- 
vestigate the cases / = 0,p = 0.65 and / = 2/5, p = 0.65. 
As shown in Figs. 7 (b) and (c), the critical current I c 
and critical exponent 1/8 for / = Q,p = 0.65 are deter- 
mined to be 0.03875 ± 0.0005, 2.24 ± 0.02 respectively, 
for / = 2/5, p = 0.65, the result is I c = 0.034 ± 0.001, 
1/8 = 2.29 ±0.01. 

We then draw the scaling plots according to Eq. ([5]). 
By one-parameter tuning (only (3), we get the best col- 
lapses of data in the regime I < I c with (3 = 2.61 ± 0.02 
and 2.28 ± 0.02 for / = 0,p = 0.86 and / = 0,p = 0.65 
respectively, which are shown in Figs. 8(a) and (b). For 
/ = Q,p = 0.86, this curve can be fitted by S(x) = 
0.0994exp(1.9a;), combined with the relation (38 — 1.55, 
suggesting a non-Arrhenius creep motion. However, for 
the strongly site-diluted system with f = 0,p = 0.65, 
the scaling curve can be fitted by S(x) = 0.037exp(0.5x), 
combined with the relation (38 « 1.0, indicative of an 
Arrhenius creep motion. Interestingly, as displayed in 
Fig. 8(c) for / = 2/5, p = 0.65, the exponent (3 is 
found to be 2.30 ± 0.02, which yields (38 w 1.0. The 
scaling curve in the regime I < I c can be fitted by 
S(x) = 0.105exp(0.25a;). These two combined facts sug- 
gest an Arrhenius creep motion in this case. 

It is worthwhile to note that both the FTPT and the 
LTCM for strongly disordered JJA's (p — 0.65) with and 
without frustration are very similar. The I-V curves 
in low temperature for all three cases can be described 
by V oc T 1 / 5 expL4(l - I c /I)/T^ s ], this is one of the 
main characteristics of glass phases [Tt], [Hj], while the 
I-V traces for KT-type phases can be fitted to V oc I a . 
Hence, by the scaling ansatze in Eq. ([5]), we have pro- 
vided another evidence for the existence of non-KT-type 
phases in the low-temperature regime for these three 
cases (/ = 0,p = 0.86;/ = Q,p = 0.65;/ = 2/5, p = 
0.65). 
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FIG. 8: (a) Scaling plot for / = 0,p = 0.86 with I c = 0.302, 
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IV. SUMMARY 

To explore the properties of various critical phenom- 
ena in site-diluted JJA's, we performed large-scale simu- 



lations for two typical percolative strengths p = 0.86 and 
p = 0.65 as in a recent experimental work 0]. We in- 
vestigated the FTPT, the ZTDT and the LTCM in these 
systems. The RSJ dynamics was applied in our work, 
from which we measured the I-V characteristics at dif- 
ferent temperatures. 

The results obtained in this work about FTPT are 
well consistent with the recent experimental findings in 
Ref. 0] and are inconsistent with the earlier experimen- 
tal study in Ref. [l[ , possibly due to the large noise in the 
measurement of voltage in Ref. [l[ (larger than 0.2nv), 
which was considerably reduced in the experiments by 
Yun et al 0] . The evidence for non-KT-type phase tran- 
sition was revealed by two different scaling ansatzes (Eq. 
(O and Eq. ©). Our results also shed some light on the 
various phases and the phase transitions where the differ- 
ent divergent correlations at various disorder strengths 
were suggested, and the critical exponents were evalu- 
ated. Furthermore, the results in this paper are useful 
for understanding not only the site-diluted systems, but 
also the whole class of disordered JJA's, for instance, the 
combination of two different phase transitions may exist 
in other disordered JJA's systems. 

In addition, the ZTDT and the LTCM were also 
touched. It was demonstrated by the scaling analysis 
that the creep law for / = Q,p = 0.86 is non-Arrhcnius 
type while those for / = 0,p = 0.65 and / = 2/5, p = 0.65 
belong to the Arrhenius type. It is interesting to note 
that the non- Arrhenius type creep law for / = 0,p = 0.86 
is similar to that in three-dimensional flux lines with a 
weak collective pinning [l7| • The product of the two ex- 
ponents 1.55 is also very close to 3/2 determined in Ref. 
|l7|]- For / = 0,p = 0.65 and / = 2/5, p = 0.65, the ob- 
served Arrhenius type creep law is also similar to that in 
the glass states of flux lines with a strong collective pin- 
ning as in Ref. [i"^ . Future experimental work is needed 
to clarify this observation. 
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We numerically investigate the intriguing effects produced by random percolative disorder in 
two-dimensional Josephson-junction arrays. By dynamic scaling analysis, we evaluate critical tem- 
peratures and critical exponents with high accuracy. It is observed that, with the introduction 
of site-diluted disorder, the Kosterlitz-Thouless phase transition is eliminated and evolves into a 
continuous transition with power-law divergent correlation length. Moreover, genuine depinning 
transition and creep motion are studied, evidence for distinct creep motion types is provided. Our 
results not only are in good agreement with the recent experimental findings, but also shed some 
light on the relevant phase transitions. 
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I. INTRODUCTION 

Understanding the critical behavior of Josephson- 
junction arrays (JJA's) with various disorders is always 
a challenging issue and has been intensely studied in re- 
cent years [l|-[I(|. However, the properties of different 
phases and various phase transitions are not well under- 
stood. Josephson-junction arrays gives an excellent re- 
alization to both two-dimensional (2D) XY model and 
granular High-T c superconductors [11] . As we know, the 
pure JJA's undergoes the celebrated Kosterlitz-Thouless 
(KT) phase transition from the superconducting state to 
the normal one, this transition is driven by the unbinding 
of thermally activated topological defects [II • When the 
disorder is introduced, the interplays among the repulsive 
vortex-vortex interaction, the periodic pinning potential 
caused by the discreteness of the arrays, and the defects 
produced by the disorder provide a rich physical picture. 

In site-diluted JJA's, the crosses around the randomly 
selected sites are removed from the square lattice. Since 
it is a representative model for realizing the irregular 
JJA's systems, how the percolation influences the physi- 
cal properties of JJA's has attracted considerable atten- 
tion [H, [3, [H 0, Q ■ Harris et al introduced random per- 
colative disorder into Nb-Au-Nb proximity-coupled junc- 
tions, the current- voltage (I-V) characteristics were mea- 
sured and the results demonstrated that the only dif- 
ference of the phase transition compared with that in 
ideal JJA's system is the decrease of critical tempera- 
ture, while the transition type still belongs to the KT 
one with the disorder strength spanning from p = 0.7 
to p = 1.0 (here 1 — p is the fraction of diluted sites) 
However, in a recent experiment, Yun et al showed 
that the phase transition changes into a non-KT-type one 
when the disorder strength increases to a moderate value 
(p = 0.86) [9(. Therefore, the existence of the KT-type 
phase transition in site-diluted JJA's remains a topic of 
controversy, the nature of these phase transitions and the 
various phases is not clear. 

On the other hand, much effort has been devoted to 



the zero-temperature depinning transition (ZTDT) and 
the related low-temperature creep motion (LTCM) both 
theoretically [II [3 [II and numerically [jj, [13, El in a 
large variety of physical problems, such as charge density 
waves [II , random- field Ising model fill ] , and flux lines in 
type-II superconductors [13, El- Since the non- linear dy- 
namic response is a striking problem, there is increasing 
interest in its properties and characteristics, especially 
in the flux lines of type-II superconductors [lT], El • In 
a recent numerical study on the three-dimensional glass 
states of flux lines, Arrhcnius creep motion was observed 
at a strong collective pinning, while the non-Arrhenius 
creep motion was demonstrated at a weak collective pin- 
ning [I3|. 

In this work, we numerically study the finite- 
temperature phase transition (FTPT) in site-diluted 
JJA's at different percolative disorder strengths, the 
ZTDT and the LTCM are also investigated. The outline 
of this paper is as follows. Sec. |TT] describes the model 
and the numerical method briefly. In Sec. IIII1 we present 
the main results, where some discussions are also made. 
Sec. IIVI gives a short summary of the main conclusions. 



II. MODEL AND SIMULATION METHOD 

JJA's can be described by the 2D XY model on a 
simple square lattice, the Hamiltonian reads 0, [20| 



H = - 



E 

<i,j> 



Jij cos( 



<t>3 ~ A a), 



(1) 



where the sum is over all nearest neighboring pairs on 
a 2D square lattice, Jy denotes the strength of Joseph- 
son coupling between site i and site j, fa specifies the 
phase of the superconducting order parameter on site i, 
and Aij — (27r/$o) / A • dl is the integral of magnetic 
vector potential from site i to site j , <&o denotes the flux 
quantum. The direct sum of around an elementary 
plaquette is 27r/, with / the magnetic flux penetrating 
each plaquette produced by the uniformly applied field, 



2 



which is measured in unit of <I>o- / = and / = 2/5 are 
the focuses of this paper. The system sizes are selected 
as 128 x 128 for / = and 100 x 100 for / = 2/5, where 
the finite size effects are negligible. We introduce the 
site-diluted disorder similar to the previous experiments 
P, Q . We first select the diluted sites randomly with the 
probability 1 — p, then remove the nearest four bonds 
around the selected sites from the lattice. The distribu- 
tions of the diluted sites are the same for all the samples 
considered. The percolative threshold concentration p c 
is about 0.592 (Hf. 

The resistivity-shunted-junction (RSJ) dynamics is in- 
corporated in the simulations, which can be described as 

turn 

<t h ^-^ ■ ■ dH ^-^ 

— 2J<& - <f>j) = -qT- + Jex,i - }^ na, (2) 
j 1 j 

where a is the normal conductivity, J ex ^ refers to the 
external current, r/ij denotes the thermal noise current 
with < T)ij(t) >= and < rj i3 (t)rj i: j(t r ) >= 2ak B T6(t - 
t). 

The fluctuating twist boundary condition is applied in 
the xy plane to maintain the current, thus the new phase 
angle 9i = <fii + ri ■ A (A = (A x , A y ) is the twist variable) 
is periodic in each direction. In this way, supercurrent 
between site i and site j is given by J'^j = Jijsin(9i — 
8j — Aij — rij • A) , and the dynamics of A Q can be written 
as 

<i,j>a 

where a denotes the x or y direction, the voltage drop 
in a direction is V = — LA a . For convenience, units are 
taken as 2e = h = Jq = a = ks = 1 in the following. 
Above equations can be solved efficiently by a pseudo- 
spectral algorithm due to the periodicity of phase in all 
directions. The time stepping is done using a second- 
order Runge-Kutta scheme with At = 0.05. Our runs are 
typically (4 — 8) x 10 7 time steps and the latter half time 
steps are for the measurements. The detailed procedure 
in the simulations was described in Ref . [2(J, H2| . In this 
work, a uniform external current I along x direction is 
fed into the system. 

Since RSJ simulations with direct numerical integra- 
tions of stochastic equations of motion are very time- 
consuming, it is practically difficult to perform any se- 
rious disorder averaging in the present rather large sys- 
tems. Our results are based on one realization of disor- 
der. For these very large samples, it is expected to exist 
a good self-averaging effect, which is confirmed by two 
additional simulations with different realizations of dis- 
order. This point is also supported by a recent study 
of JJA's by Um et al Q. In addition, simulations with 
different initial states are performed and the results are 
nearly the same. Actually, the hysteretic phenomenon is 



usually negligible in previous RSJ dynamical simulations 
on JJA's [3, Q ■ For these reasons, the results from simu- 
lations with a unique initial state (random phases in this 
work) are accurate and convincing. 

III. RESULTS AND DISCUSSIONS 
A. Finite temperature phase transition 

The I-V characteristics are measured at different dis- 
order strengths and temperatures. At each temperature, 
we try to probe the system at a current as low as possi- 
ble. To check the method used in this work, we investi- 
gate the I-V characteristics for / = 0,p = 1.0. As shown 
in Fig. 1(a), the slope of the I-V curve in log-log plot 
at the transition temperature T c (« 0.894) is equal to 
3, demonstrating that the I-V index jumps from 3 to 1, 
consistent with the well-known fact that the pure JJA's 
experiences a KT-type phase transition at T c w 0.894. 
Figs. 1(b) and (c) show the I-V traces at different dis- 
order strengths in unfrustrated JJA's, while Fig. 1(d) is 
for / = 2/5, p = 0.65. It is clear that, at lower temper- 
atures, R — V/I tends to zero as the current decreases, 
which follows that there is a true superconducting phase 
with zero linear resistivity. 

It is crucial to use a powerful scaling method to an- 
alyze the I-V characteristics. In this paper, we adopt 
the Fisher-Fisher-Huse (FFH) dynamic scaling method, 
which provides an excellent approach to analyze the 
superconducting phase transition (23|. If the prop- 
erly scaled I-V curves collapse onto two scaling curves 
above and below the transition temperature, a continu- 
ous superconducting phase transition is ensured. Such a 
method is widely used @, [24[ , the scaling form of which 
in 2D is 

v = ic z MiO, (4) 

where ip + r_\(x) is the scaling function above (below) T c , 
z is the dynamic exponent, £ is the correlation length, 
and V - I z+1 at T = T c . 

Assuming that the transition is continuous and char- 
acterized by the divergence of the characteristic length 
f ~ |T - T c \~ v and time scale t ~ £ z , FFH dynamic 
scaling takes the following form 

(V/I)\T - T c \- Z » = i, ± (I\T - T c \~»). (5) 

On the other hand, to certify a KT-type phase tran- 
sition in JJA's, a new scaling form [25j is proposed as 
follows 

(I/T)(I/V)V*=P ± {IZ/T), (6) 

which can be derived directly from Eq. ^ after some 
simple algebra. The correlation length of KT-type phase 
transition above T c is well defined as £ ~ e ( c /l T - T d) 1/2 
and Eq. ^ is rewritten as 

(I/T^I/V) 1 ^ = P+(Ie^ T - T ^ 1/2 /T). (7) 
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FIG. 1: I-V characteristics for different frustrations and dis- 
order strengths. The dash lines are drawn to show where 
the phase transition occurs, the slopes of which are equal to 
z + 1, z is the dynamic exponent. The transition temperature 
and dynamic exponent for (a) are well consistent with the 
well-known result, i. e., T c = 0.894, z = 2.0, for (b),(c),(d) 
are well consistent with those determined by FFH dynamic 
scaling analysis. Solid lines are just guide to eyes. 



We perform the dynamic scaling analysis at a strong 
disorder (p — 0.65) in unfrustrated system (/ = 0). Us- 
ing T c = 0.24±0.01, z = 1.2±0.02 and v = 1.0±0.02, an 
excellent collapse is achieved according to Eq. ([5]) , which 
is shown in Fig. 2. In addition, all the low-temperature 
I-V curves can be fitted to V ~ iexp(— with 
/i = 0.9 ~ 1.1. These results certify a continuous super- 
conducting phase with long-range phase coherence. The 
critical temperature for such a strongly disordered system 
is very close to that in 2D gauge glass model (T c = 0.22) 

For / = 0,p = 0.86, we first still adopt the scaling 
form in Eq. ([5]) to investigate the I-V characteristics. 
As displayed in Fig. 3, we get a good collapse for T <T C 
with T c = 0.58 ± 0.01, z = 2.0 ± 0.01 and u = lA± 0.02, 
demonstrating a superconducting phase with power-law 
divergent correlation for T < T c . Note that the collapse 
is poor for T > T c , implying that the phase transition 
is not a completely non-KT-type one. Next, we use the 
scaling form in Eq. to analyze the I-V data above 
T c . Interestingly, using T c = 0.58 and z = 2.0 determined 
above, a good collapse for T > T c is achieved, which is 
shown in Fig. 4. That is to say, the I-V characteristics 
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FIG. 2: Dynamic scaling of I-V data at various temperatures 
according to Eq. ([5]) for / = 0,p — 0.65. 



at T < T c are like those of a continuous phase transition 
with power-law divergent correlation length while at T > 
T c are like those of KT-type phase transition, which are 
well consistent with the recent experimental observations 

1- 

To make a comprehensive comparison with the experi- 
mental findings in Ref . 0] , we also investigate the FTPT 
in frustrated JJA's (/ = 2/5) at a strong site-diluted dis- 
order (p = 0.65). As shown in Fig. 5, a superconducting 
phase transition with power-law divergent correlation is 
clearly observed. 

As is well known, non-KT-type phase transition in 
frustrated systems is a natural result. However, it is in- 
triguing to see that in unfrustrated systems, one may 
ask what our results really imply and what is the mecha- 
nism. It has been revealed that in the presence of a strong 
random pinning which is produced by random site dilu- 
tions, a breaking of ergodicity due to large energy barrier 
against vortex motion may allow enough vortices to ex- 
perience a non-KT-type continuous transition (28j . 

Interestingly, we recover the phenomena in experi- 
ments by the present model and give some insights into 
the FTPT. More information on the low-temperature 
phase calls for further equilibrium Monte Carlo simu- 
lations as in Refs. [27| . Table I summarizes the crit- 
ical temperatures at different frustrations and disorder 
strengths. One can find that, the critical temperature in 
unfrustrated system decreases with increasing number of 
diluted sites. 



TABLE I: Summary of T c . 





f=0 


£=2/5 


p=0.95 


0.85(2) 


0.16(2) 


p=0.86 


0.58(1) 


0.13(1) 


p=0.7 


0.27(2) 


0.12(1) 


p=0.65 


0.24(1) 


0.14(1) 



The systems considered in our work are site-diluted 
JJA's, which are not the same as bond-diluted JJA's in 
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FIG. 3: Dynamic scaling of I-V data at various temperatures 
according to Eq. © for / = 0,p = 0.86, T < T c . Solid lines 
are just guide to eyes. 



FIG. 5: Dynamic scaling of I-V data at various temperatures 
according to Eq. ([5]) for / = 2/5, p = 0.65. 
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FIG. 4: Dynamic scaling of I-V data at various temperatures 
according to Eq. © for / = 0,p = 0.86, T > T c . Solid lines 
are just guide to eyes. 



Ref. H |H . In bond-diluted systems the diluted bonds 
are randomly removed, while in the site-diluted systems, 
the diluted sites are randomly selected, then the nearest 
four bonds around the selected sites are removed. Al- 
though the JJA's in Ref. [1, H[ and the present work are 
diluted in different ways, it is interesting to note that the 
obtained exponents in FTPT are very close, possibly due 
to the similar disorder effect produced. 



B. Depinning transition and creep motion 

Next, we turn to the ZTDT and the LTCM for the 
typical site-diluted JJA's systems mentioned above. De- 
pinning can be described as a critical phenomenon with 
scaling law V ~ (7 — ic) , demonstrating a transition 
from a pinned state below critical driving force I c to 
a sliding state above I c . The (I — I c ).vs.V traces at 
T = for / = 0,p = 0.86; / = 0,p = 0.65 and 
/ = 2/5, p = 0.65 are displayed in Fig. 6, linear-fittings 
of Log(I — I c ).vs.LogV curves are also shown as solid 
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FIG. 6: (a)/y characteristics for / = 0,p = 0.86 with I c = 
0.302 ± 0.005, P = 2.62 ± 0.1. (bjlV characteristics for / = 
0,p = 0.65 with I c = 0.039 ± 0.001, = 2.37 ± 0.1. (c)IV 
characteristics for / = 2/5, p = 0.65 with I c = 0.035 ± 0.002, 
f3 = 2.27 ±0.05. 
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satisfy the following scaling ansatz [17j, [26[ 

V(T, I) = T^ S S[{1 - IJI)^ 1 / 135 ]. 



(8) 
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FIG. 7: (a)LogT-LogV curves for / = 0,p = 0.86 around I c 
with I c = 0.302 ± 0.001, 1/8 = 1.688 ± 0.001. (b)LogT-LogV 
curves for / = 0,p = 0.65 around I c with I c = 0.03875 ± 
0.0005, 1/5 = 2.24 ±0.02. (c)LogT-LogV curves for / = 
2/5, p = 0.65 around I c with I c = 0.034 ±0.001, 1/8 = 2.29 ± 
0.01. 



lines. As for f = 0,p = 0.86, the depinning exponent (3 
is determined to be 2.62 ± 0.1 and the critical current I c 
is 0.302 ± 0.005, while for the cases f = 0,p= 0.65 and 
/ = 2/5, p = 0.65, the depinning exponents are evaluated 
to be 2.37 ±0.1 and 2.27 ±0.05 with the critical currents 
I c = 0.039 ± 0.001 and I c = 0.035 ± 0.002, respectively. 

When the temperature increases slightly, creep mo- 
tions can be observed. In the low-temperature regime, 
the I-V traces are rounded near the zero-temperature 
critical current due to thermal fluctuations. Fisher first 
suggested to map such a phenomenon for the ferromag- 
net in magnetic field where the second-order phase tran- 
sition occurs [29( | . This mapping was then extended to 
the random-held Ising model [161 ] and the flux lines in 
type-II superconductors [llj • For the flux lines in type-II 
superconductors, if the voltage is identified as the order 
parameter, the current and the temperature are taken as 
the inverse temperature and the field respectively, anal- 
ogous to the second-order phase transition in the ferro- 
magnet, the voltage, current and the temperature will 



Where S(x) is a scaling function. The relation V(T,I = 
I c ) = S f (0)r 1//(5 can be easily derived at I — I c , by which 
the critical current I c and the critical exponent 8 can be 
determined through the linear fitting of the LogT-Logy 
curve at I c . 

The LogT-LogV curves are plotted in Fig. 7(a) for 
/ = 0,p = 0.86. We can observe that the critical current 
is between 0.3 and 0.32. In order to locate the criti- 
cal current precisely, we calculate other values of voltage 
at current within (0.3,0.32) with a current step 0.01 by 
quadratic interpolation [26j]. Deviation of the T-V curves 
from the power law is calculated as the square deviations 
SD = ^2[V(T) — y(T)] 2 between the temperature range 
we calculated, here the functions y(T) = C1T~ C2 are 
obtained by linear fitting of the LogT-LogV/ curves. The 
current at which the SD is minimum is defined as the 
critical current. The critical current is then determined 
to be 0.302 ± 0.001. Simultaneously, we obtain the expo- 
nent 1/8 = 1.688 ± 0.001 from the slope of LogT-LogV 
curve at I c = 0.302. The similar method is applied to in- 
vestigate the cases / = 0,p = 0.65 and / = 2/5, p = 0.65. 
As shown in Figs. 7 (b) and (c), the critical current I c 
and critical exponent 1/8 for / = Q,p = 0.65 are deter- 
mined to be 0.03875 ± 0.0005, 2.24 ± 0.02 respectively, 
for / = 2/5, p = 0.65, the result is I c = 0.034 ± 0.001, 
1/8 = 2.29 ±0.01. 

We then draw the scaling plots according to Eq. ([5]). 
By one-parameter tuning (only (3), we get the best col- 
lapses of data in the regime I < I c with (3 = 2.61 ± 0.02 
and 2.28 ± 0.02 for / = 0,p = 0.86 and / = 0,p = 0.65 
respectively, which are shown in Figs. 8(a) and (b). For 
/ = Q,p = 0.86, this curve can be fitted by S(x) = 
0.0994exp(1.9a;), combined with the relation (38 — 1.55, 
suggesting a non-Arrhenius creep motion. However, for 
the strongly site-diluted system with f = 0,p = 0.65, 
the scaling curve can be fitted by S(x) = 0.037exp(0.5x), 
combined with the relation (38 « 1.0, indicative of an 
Arrhenius creep motion. Interestingly, as displayed in 
Fig. 8(c) for / = 2/5, p = 0.65, the exponent (3 is 
found to be 2.30 ± 0.02, which yields (38 w 1.0. The 
scaling curve in the regime I < I c can be fitted by 
S(x) = 0.105exp(0.25a;). These two combined facts sug- 
gest an Arrhenius creep motion in this case. 

It is worthwhile to note that both the FTPT and the 
LTCM for strongly disordered JJA's (p — 0.65) with and 
without frustration are very similar. The I-V curves 
in low temperature for all three cases can be described 
by V oc T 1 / 5 expL4(l - I c /I)/T^ s ], this is one of the 
main characteristics of glass phases [Tt], [Hj], while the 
I-V traces for KT-type phases can be fitted to V oc I a . 
Hence, by the scaling ansatze in Eq. ([5]), we have pro- 
vided another evidence for the existence of non-KT-type 
phases in the low-temperature regime for these three 
cases (/ = 0,p = 0.86;/ = Q,p = 0.65;/ = 2/5, p = 
0.65). 
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FIG. 8: (a) Scaling plot for / = 0,p = 0.86 with I c = 0.302, 
1/5 = 1.688 and (35 = 1.55. (b) Scaling plot for / = 0,p = 
0.65 with Ic = 0.03875, 1/5 = 2.24 and (35 » 1.0. (c)Scaling 
plot for / = 2/5, p = 0.65 with I c = 0.034, 1/5 = 2.29 and 
/3<5 « 1.0. 



IV. SUMMARY 

To explore the properties of various critical phenom- 
ena in site-diluted JJA's, we performed large-scale simu- 



lations for two typical percolative strengths p = 0.86 and 
p = 0.65 as in a recent experimental work 0]. We in- 
vestigated the FTPT, the ZTDT and the LTCM in these 
systems. The RSJ dynamics was applied in our work, 
from which we measured the I-V characteristics at dif- 
ferent temperatures. 

The results obtained in this work about FTPT are 
well consistent with the recent experimental findings in 
Ref. 0] and are inconsistent with the earlier experimen- 
tal study in Ref. [l[ , possibly due to the large noise in the 
measurement of voltage in Ref. [l[ (larger than 0.2nv), 
which was considerably reduced in the experiments by 
Yun et al 0] . The evidence for non-KT-type phase tran- 
sition was revealed by two different scaling ansatzes (Eq. 
(O and Eq. ©). Our results also shed some light on the 
various phases and the phase transitions where the differ- 
ent divergent correlations at various disorder strengths 
were suggested, and the critical exponents were evalu- 
ated. Furthermore, the results in this paper are useful 
for understanding not only the site-diluted systems, but 
also the whole class of disordered JJA's, for instance, the 
combination of two different phase transitions may exist 
in other disordered JJA's systems. 

In addition, the ZTDT and the LTCM were also 
touched. It was demonstrated by the scaling analysis 
that the creep law for / = Q,p = 0.86 is non-Arrhcnius 
type while those for / = 0,p = 0.65 and / = 2/5, p = 0.65 
belong to the Arrhenius type. It is interesting to note 
that the non- Arrhenius type creep law for / = 0,p = 0.86 
is similar to that in three-dimensional flux lines with a 
weak collective pinning [l7| • The product of the two ex- 
ponents 1.55 is also very close to 3/2 determined in Ref. 
|l7|]- For / = 0,p = 0.65 and / = 2/5, p = 0.65, the ob- 
served Arrhenius type creep law is also similar to that in 
the glass states of flux lines with a strong collective pin- 
ning as in Ref. [i"^ . Future experimental work is needed 
to clarify this observation. 
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